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Abstract 

Validity of the mean-field approach to open system dynamics in the optical cavity 
system is examined. It is rigorously shown that the mean-held approach is justified in 
the thermodynamic limit. The result is applicable to nonequilibrium situations, e.g. 
the thermal reservoirs may have different temperatures, and the system may be subject 
to a time-dependent external field. The result of this work will lead to further studies 
on macroscopic open quantum systems. 

1 Introduction 

Quantum dynamics in many body systems is a very important topic for many branches 
of physics. In particular, recent studies in the field of statistical physics have focused 
on dynamical problems such as thermalization [1-5], phase transitions induced by the 
parameter quench [6-9], and nonequilibrium phase transitions [10-12]. Theoretically, it is 
difficult to precisely describe the time evolution of a many body system. Although recent 
works succeeded to treat the exact time evolution and observe relaxation processes in some 
integrable systems [13-15], we should be content with some approximate treatments in 
general. 

One of the important approximate approaches is the MF approximation for not only 
equilibrium statics but also nonequilibrium dynamics. In this approximation, the iV-body 
state vector or the A^-body density matrix is approximated by the product state. Interest- 
ingly, the MF approach predicts some remarkable results such as the absence of thermal- 
ization and purely dynamical phase transitions in some isolated quantum systems [9, 10]. 
Also in open quantum systems, nonequilibrium phase transitions have been studied by the 
MF approach [16-18], and see also Ref. [19] for the recent result. 

Some of these predictions might be the peculiarity of the MF theory and not universal 
in general. But, importantly, there are several models in which the MF treatment becomes 
exact in some ideal limit [20-23] . Thus the MF theory is completely reliable as long as such 
an ideal limit is considered to be realized. For example, in a quantum spin system with 
a global coupling (the infinite-range interaction), it is known that the quantum dynamics 
is exactly described by the MF dynamical equation (or Hartree equation) [20]. Another 
example is the iV-body bosonic system interacting via the two-body potential of the scat- 
tering length a, in which the time-dependent Gross-Pitaevskii equation, which is regarded 
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as the MF theory of the dynamics of the Bose-Einstein condensate, becomes exact in the 
limit of N — >■ oo with a fixed value of Na [21]. Because of the existence of such ideal limits, 
the MF approximation takes a special position among other approximate approaches. 

Therefore, we should not regard the MF approach as a mere crude approximation. 
Rather it is important to clarify the condition under which the MF approach is justified 
and extend the possibility of experimental realization of remarkable predictions of the MF 
approach. 

This article is mainly devoted to investigating the quantum dynamics of the system 
consisting of TV two-level atoms interacting with a single quantized boson mode in an 
optical cavity. This system is described by the Dicke model [24]. We consider its extension, 
that is, we allow each element of the system, i.e. single cavity mode of photons and iV 
two-level atoms, to be attached to its own environment. Each environment is a Hamilton 
system, and may or may not be a large system. If one environment has a continuous energy 
spectrum, which means that the environment is an infinitely large system, it plays the role 
of a dissipative thermal bath. 

The MF approach of this model predicts interesting dynamical or nonequilibrium phase 
transitions. In the isolated Dicke model (there is no environment), it was argued that 
dynamical phase transitions occur by the MF approach [9]. In open systems (the Dicke 
model interacting with dissipative environments), under the time-periodic driving field, the 
MF dynamics combined with the method of the Born-Markov quantum master equation [25] 
predicts some nonequilibrium phase transitions including optical bistabilities [17, 18] which 
were observed in experiments [26]. It has been believed that the MF approach is justified 
in the limit of N — > oo because in this model N two-level atoms interact with a common 
cavity mode of photons, and it is expected that each atom feels only the MF produced 
by the other atoms via radiation and absorption of cavity photons. In this situation, the 
correlation between two atoms or between an atom and cavity photons will be not so 
important. However there has not been rigorous proof of this expectation yet. The aim of 
this work is to provide its rigorous proof. 

The organization of this paper is as follows. In Sec. 2, we introduce the proof of the 
fact that the MF theory is justified in infinite-range spin models. This result is already 
known [20], but we review it because it is the simplest case where the MF approach is 
justified rigorously. In Sec. 3, the model considered in the present work is explained. In 
Sec. 4, some preliminaries are made. We introduce the two notions important later, i.e. the 
Q-representation of the density matrix and restriction of the class of physical observables. 
In Sec. 5, the main result of this paper is stated, and some comments on this main result 
are mentioned. In Sec. 6, the main result is proven. 

2 Exactness of the mean-field dynamics in infinite-range spin 
models 

Before studying the optical cavity system, let us analyze the simplest case where the MF 
dynamics becomes exact in the thermodynamic limit. We consider the N spin system 
with infinitely long-range interactions. The Hilbert space for ith spin is denoted by Hi- 
Let hi = h a bounded self-adjoint operator acting on Hi, and let Vij = V (i ^ j) also a 
bounded self-adjoint operator acting on % QHj. We put Vu = 0. The Hamiltonian of the 
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spin system with infinite-range two-body interactions is generally written as 



1 N 



(1) 



(2) 



(3) 



The operator A is said to be bounded if ||-A||oo < +°°- 

The iV-spin density matrix at time t is denoted by pN,t, and its /c-marginal reduced 



We introduce the p-norm of an operator A acting on ® i=1 "Hi as 

Pllf := (Tr 1 , 2 ,..., k \A\ P ) 1/P ■ 
In particular, the operator norm is defined as 

\\A\\oo ■= sup ■ 



density matrix by 7]^ which is defined as 

7jv,t := lr fc+i,fc+2,...,AfPAf,t- 
The initial condition is assumed to satisfy 

VfceN, 



lim 

7V->oo 



?Ar,0 To 



(fc) 
1 



(4) 



(5) 



that is, reduced density matrices are initially factorized. Here 70 is a single spin density 
matrix. 

The MF theory assumes that reduced density matrices are always factorized, 



T^i ~ Tmf,* = 7t®* VfcGN. 
The single spin density matrix 74 obeys the Hartree equation, 



(6) 



d 



-idp-yt ~ *Th4?7t ® Tt 
-«[/ii,7t] +Tr 2 [Fi2,7* ®7f]- 



(7) 



The Liouville operators are defined as £^(-) := [/tj, (•)] and := [Vij, (•)]. 

The statement of exactness of the MF dynamics is that if the initial state satisfies 
Eq. (5), 



lim 

N-xx 



(*) _ 
!N,t It 



(k) 



= 



(8) 



for any fixed k € N and i > 0. 1 

From Eq. (8) it is immediately verified, by using the inequality ||AB||^ < 1 1 ^4 1 1 oo 1 1 [ | ^ fc ^ , 
that 



lim 



N,t 



MF.t 



0. 



(9) 



1 Exactness of the MF dynamics does not mean p N>t — 7®^. Actually, it does not hold even in infinite- 
range spin models. 
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for any bounded operator acting on ^)^ =1 Hi- Here ( ■ ) jv,t := Tr(-)pjv,t ano - (")mf, 
Tr 1)2 ,..., fc (.) 7 Si, t . 

We follow Ref. [22] for the proof of Eq. (8). We start with the Liouville equation 
d 



d,t 



PN,t 



/ N N \ 

\»=i u / 



PN,t- 



(10) 



By tracing out over 0^ fc+1 we obtain the hierarchical equations (BBGKY hierarchy 2 ) 



& hffj} [22], 

i=i 

The super-operators and w( fc ) are defined by 



The MF solution 7^ t = ^/f k satisfies the equation 



(11) 



(12) 



(13) 



which is obtained by formally taking the limit of N — >■ 00 in Eq. (11). 

It should be remarked that only the fact that Eq. (13) is obtained from Eq. (11) by 
formally taking the limit of N — > 00 does not ensure that the MF dynamics is exact in 
this limit. Because Eq. (11) forms a coupled chain of N equations of motion, the terms 
proportional to 1/iV in Eq. (11) might be amplified in the chain and have nonnegligible 
contribution. In particular, for a very large k < N, the last two terms in RHS of Eq. (11) 
is no longer small. We must prove that these terms actually do not trouble. 

The next step is to construct the Duhamel series [22] by formally integrating Eq. (11). 
We define 



U^ k) := exp 



i=i 



Then we obtain 



J 



(14) 
(15) 



"Bogoliubov-Born-Green-Kirkwood-Yvon hierarchy 
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We substitute this expression iteratively to the terms in RHS not proportional to 1/N (the 
second term of RHS in the above equation), and repeat this procedure L times, then we 
obtain the following so called the Duhamel series expansion, 

l=Q J J ■/ 

+ fdtr Pdt 2 ... /^W^*^ 
Jo Jo Jo 

n (k+l-l) v (k+l-l) (k+l-l) 



+ N 



i=i 



r 7 (*+i-l) W (fc+J-l)„(*+0 



(16) 



The Duhamel series of the MF density matrix is written 

A = E /" *i ■ ■ ■ r _1 ^A^K^ 1 ' • • • w^- 1 



u t l 7mf,o 



+ *i jf ■ ■ ■ jf " ^W-l ^'^(W) . . . U^W^ L 'V£il (17) 



From these expressions, we shall show ||7^ — TMFtlli^ — ^ in the limit of N 
From Eqs. (16) and (17), we obtain 



(k) _ (fc) 



(fc) 

< A! + A 2 + A 3 + At, 



oo. 



(18) 



where A; L (i = 1, 2, 3,4) are given by 



Ai = E f dti P dt 2 ... P- 1 dt t U^W^U^W^ 
i=o ^° ^° ^° 



• 1 u t, \ in,o Tmf,o 



(fe) 



A, 



jo Jo Jo 

r/ (fc+L-l) w ( fe+L _i) / (fe+L) (fc+m 
• • • U t L -i-t L VV \'N,t L IMF,t L ) 



(fc) 
1 
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and 



A3 



I L ft rh rh-i 



dtl 



tl-*2 



r7 (fc+J-l) v (fc+Z-l) 



(k) 
i 



A 4 = V - dh dt 2 ... dti 

j^i ^ Jo Jo Jo 



r/ (*+i-i) w (fc+i-i) (fc+i) 



(fc) 
i 



Now we evaluate Aj. Since J7 t ^ is a unitary operator, it does not change the norm, 
||^( fe )(.)||( fe ) = ||(.)||J fe) . We further use the inequality \\AB\\[ k) < PHo^BH^. As a result, 
Ai is evaluated, by recalling the definitions of and W^ k \ as 



1=0 



i( k + l-l)\ 
(fc-1)! 



(fc+0 (fc+o 
7at,o 7mf,o 



(fc+0 
l 



By the assumption for the initial condition, Eq. (5), ii^OasiV^ooata fixed L. 
Similarly, the upper bound of A2 becomes 

t 



A 2 <-{2\\V\\ 00 



(fc + L-1)! 



< 



(fc + L-1)! 
L!(fc-1)! 



(fc-1)! 
2||^||oot) L 



(fc+j) _ (fc+J) 



(fc+0 
1 



In the last inequality, we used 
2 k+L -\ we have 



(k+i) _ (k+l) 
lN,t L 7mf,* l 



{k+l) 



< 2. Since (k + L - l)l/[L\(k - 1)!] < 



a 2 < 2 fc (4||y|| 00 t) L . 



If we restrict the time to t < to := 1/ (8|| V||oo) (this restriction will be removed later), we 
have A 2 < 2 h ~ L . If we take the limit of L — > 00 after N — > 00, both Ai and A 2 converge to 
zero. 

Similarly, we can show that A3 and A4 also converge to zero in the same limit. The 
upper bound of A3 is evaluated as 



= ±E(2\\v\uy(k + i-i)( k + l - l)l 



N ^ v ""' '" M " / v ~ ' " U(k-1)\ 



(19) 



1=1 



If we take the limit of N — > 00 first, obviously A3 goes to zero. A4 is also bounded from 
above by Eq. (19). Therefore, it has been proven that limL_>. 00 lim7v_ 5 . 00 Aj = for the time 
interval < t < to- 
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The restriction of to is not essential. Since 




satisfies the assumption of Eq. (5), we 



can start with |t]v t j a ^ time to as a new initial condition. By applying the same argument 
above, it is concluded that Eq. (8) holds in the time interval < t < 2to- Repeating this 
argument, we can extend to — > oo, and Eq. (8) is proved for an arbitrary time t > 0. 

The above proof relies on the boundedness of the operators. In optical systems, however, 
the above proof cannot be used as it is since the creation and the annihilation operators 
of bosons are unbounded. We construct the proof for such a situation in the following 
sections. 

3 Model 

We consider the Schrodinger dynamics of the generalized Dicke model given by 



H(t) = 


TV 
j=0 


'H So (t) = 


w p o t o+ 0V^(t)(a + a t ), 


H Sj = 


N 

"aJ2 S J l? = l,2,...,A0, 

i=l 


H B = 

< 


^2u r blb r , 

r 


H h = 


^ A r (6 r a^ + b\.a), 

r 


V = 

1. 


N 

V j=i 



(21) 



In (t) , the index j = corresponds to the single mode of cavity photons driven by the 
external field £(t), and j = 1,2,...,N corresponds to the ensemble of N two-level atoms. 
Cavity photons (i?s ) interact with the free Bose field Hb through the coupling H\ . 
Each two- level atom H^j, j = 1, 2, ... N, interacts with an arbitrary Hamilton system H-Q j 
through an arbitrary interaction Hamiltonian H\ j . We call the system of the Hamiltonian 
Hb 3 "j-environment" (j = 0,1,..., iV), and call the composite system described by the 
Hamiltonian Ho lj (t) + B-Q j + H\ j "j-subsystem" . We assume that H-q. + H\ j are identical 
for j = 1,2, ... ,N; the Hamiltonian H(t) is symmetric under the exchange of two indices 
i and j (i, j = 1,2, . . . ,N). 

The last term of Eq. (20) represents the interaction between the cavity photons and the 
ensemble of atoms. If we choose X,- = (S~), the Hamiltonian J2f=o {t) + V is called 
the Dicke model (Tavis-Cummings model), respectively, with a driving force £. Therefore, 
we can regard Eq. (20) as an extension of the Dicke-like model; environmental systems are 
attached to it. 

If the j-environment has a continuous spectrum, then this environment acts as a thermal 
reservoir in contact with H^^ ■ We can obtain the Hamiltonian with a continuous spectrum 
as a limiting case of the discrete spectrum. It is necessary to properly choose the coupling 
constants {A r } between cavity photons and 0-environment in order to have a well-defined 
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limit. We introduce a parameter A so that the number of eigenmodes with the frequency 
between u and uj + du is given by AD(u)dw. Thus the limit of A — > oo corresponds to 
the limit of a continuum of O-environmental modes. We assume that A r ~ X(uj r ) / y/A. 
In this case, if we consider the limit of the continuous spectrum of the O-environment, 
^ r — > A J °° dojD(cj), and therefore, ^ r |A r | 2 — )• J °° dujD(oj)\(uj) 2 =: J °° duiJ(u) does not 
depend on A. Here, J(oj) is called the spectral density [27]. We assume that 

^|A r | 2 <+oo (22) 

r 

in order not to diverge the effect of the O-environment on the dynamics of cavity photons. 
In general, for the proof of the justification of the MF treatment, A is arbitrary (it can 
even depend on N). When A <C N, the effect of the O-environment on the dynamics of the 
whole system is negligible. When A ~ N, the O-environment affects the dynamics of the 
total system. When A 3> N, the O-environment behaves as a thermal reservoir, and cause 
the dissipation. 



Mean-field dynamics 



The density matrix of the whole system, pN,t, is defined on the Hilbert space (S>iLo where 
T-Lj (j = 0, 1, 2, . . . , N) is the Hilbert space of the j-subsystem. The MF theory assumes 

(k) 

that the /c-marginal reduced density matrix x N t = Trk+i t ... t NPMF,t is in the product form, 

7S«7S^ = (^®(Pa)f fc . (23) 

(p p )t is the density matrix for the 0-subsystem. Similarly, (p a )t is a common density 
matrix for the A;-subsystem (k = 1,2, .. . , N). The density matrices (p p )t and (p a )t obey 
the following Hartree equations, 

J t (Pv)t = ~ 1 K atfl + V^Xa + « f ) + #B + #lo 

+,/Ng {[TrXMtla! + ^Ng[TrXl(p a ) t ]a} , (p p ) t ] , (24) 



j t {p,)t 



u a S{ + H Bl +H h + -jL= [[Tra(p p ) t ]xl + [Tra^p^X^ , ( Pa ) t 



(25) 



The apparent A^-dependence can be removed by the appropriate scaling in the coherent 
state representation, see Sec. 4.1. 

It is stressed that the MF theory here corresponds to neglect the correlation between 
i and j-subsystems. In general, we cannot neglect the correlation between the system of 
interest (i?s,) an d environments (Hi,). 



4 Preliminaries 

Although we cannot justify the MF theory straightforwardly in the same way as we did in 
the spin systems, we follow essentially the same course. We derive the BBGKY hierarchy 
for the reduced density matrices, and construct the Duhamel series expansion, and evaluate 
each term. However, because of unboundedness of boson operators we cannot use the trace 
norm, and the statement itself must be modified. Some preliminaries are necessary before 
presenting the main result of this work. 
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4.1 Duhamel series expansion in the coherent state representation 

Since the density matrix pN jt obeys the Liouville equation 

d 



dt 



PN,t = -i[H(t),p N>t ] =: -iCp N j, 



(26) 



we obtain the following chain of equations of motion for |7^ t |, 

k 

^(^.+iJ B ,+^i,),7g 



d (k) 

r N ; t = - % 



dt 



3=0 



N 



3 = 1 



i(N - k)Tr k+1 



(aXl +1 + atx fc+ i) , 



(fc+i) 



T N . 



(27) 



Because it is hard to see which terms are important in the limit of N — > oo, and because 
the relation to the MF theory is not obvious in this form, let us introduce the Husimi's 
"Q-representation" of the reduced density matrix [28] : 



q!S(«,{M) 



9 



JV+1 



0Va;{0V/3 r } 7 $ y^Na; {y/N/3 r } 



(28) 



Here jV is the number of modes of O-environment (j\f = ^ r 1 ~ J °° D{uj)duj). The vector 



iVa, {y^r}) denotes the coherent state on the Hilbert space T-Lq, and it satisfies 



Na,{^Nf3 r }) o , 



(29) 
(30) 



Notice that is still a matrix on the Hilbert space ®j=i H-j- 
Q$ t has the following properties: 

(i) Q { n\ > 0, 

(ii) J d 2 a (j[ J d 2 ^ Tn,2,...,kQ% = 1, 

(iii) Tr 0il ,... ifc O( fc ) 



a a' 



b r bl 



(k) 



= / ^ a (II / Tr li2i ..., fe O«(a,a*,{/3 r ,/3;})Q 



(fc) 



vr=l" 
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where is an arbitrary anti-normal ordered operator acting on ^)^ =Q %j\ namely, all 
the annihilation operators are in the left of all the creation operators. From the above 
properties, we can obtain the expectation value of any observable from Q^\- 

(k) (k) 

An advantage to use Q N t instead of t n t is that it becomes easy to see which terms are 
important and which terms are likely to be negligible in the limit of N — > oo. Indeed, the 



time evolution equation of Q N t is given by 

■i(c p (t) + 4 fc) («))oSi + ( i - w«q^ + Ivwgw , 



where 



-i£ p (t)Q^ N \ 



;— 
da 



(jj. 



(k) 

N,t 



(31) 



— i 



d 
da* 



,a + f (t) + ^2 X rPrj Q 



d 

Wr 



(co r f3 r + \ r a) Q 



(*)• 

AT.t 



(fc) 
AT,t 

d 



(k) 



(oj r p; + X* r a*) Q% 



(32) 



which corresponds to the free time evolution of the O-subsystem, and 

k k 

(fe) 



-iCf\a)Q% : 



(u a S* + H Bj + Hi.) [aXj + a*X, 

3=1 3=1 



Q 



(k) 



3 ; ' N,t 



(33) 



which represents the time evolution for k atoms and environments under the "effective 
field" [aX] + a*X^. The super-operator V (fc) and W {k) are defined by 



3=1 



3 = 1 



da* 




(34) 



(35) 



As was mentioned in Sec. 3, it is remarked that if we consider the case in which the 
O-environment acts as a thermal bath, the limit of A — > oo should be taken before N — >■ oo. 
By assumption, {A r } satisfies J2 r I'M 2 < oo in the limit of A — > oo. In this case, the state of 
the O-environment is almost unchanged, that is, (b r ) ~ \J~Nf} r ~ 0(^/N\ r ) ~ 0(yfN/K), 
if the O-environment is initially in equilibrium, so it is consistent with the interpretation 
that the O-environment is a thermal reservoir. On the other hand, as is seen in the above 
estimation, when A ~ O(N), the state of the O-environment is strongly disturbed, and it 
cannot be regarded as a thermal bath in this case. Anyway, the justification of the MF 
approach is possible for both these two cases. 

If we formally take the limit of N — > oo in Eq. (31), we obtain the following equation 
of motion: 

|qS, = -<0Cp(t) + 4 fc) (a))Qi& >4 + wWQiS?. 



(36) 
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If initially the density matrix is in the product form, 

Qmf,o = Qmf,o®(p*)o\ VfcGN 

the solution of Eq. (36) is equivalent to that of the MF dynamical equation. Therefore, we 
shall compare the solutions of Eq. (31) with those of Eq. (36). 

We define the time evolution operator due to — iC p (t) — iC^\a) as 

:= T" exp Li fdt' (£p(0 + 4 fc) («))l • (37) 

■J S 

Here T~ is the time-ordering operator (the arrow implies the direction from past to future). 

(k) 

Then we can obtain the formal integral equation for Q N t from Eq. (31) as 

Q% =u<$Q%{0) + f ' dtM> {k) Q^ ] 

J 

It is almost the same as Eq. (11). By substituting this expression into RHS iteratively L 
times, we obtain the following Duhamel series expansion: 

Q% = E f *i f r _1 ^C'^) . . . w^u^q^ 

jTq Jo Jo Jo 

r dtl r dt2 ... r _i dt L u^ k M k j> {k+ " ] ■ ■ ■ uesm^q^? 

J J 

^ •/ •/ 

- E jf * jf *» ■ ■ ■ jf"' ^^ogs'hv^') 

...^I,T 1) W<'+'-')<?S5«. (39) 
Similarly, we can obtain the Duhamel expansion for the MF density matrix Q MFt - 

oS&, t = E /* *i f ■ ■ ■ f " ^<^<+> (fc+1) . . . w^MTq ( m^ 

]Tq Jo Jo Jo 

•/ii ./ii 

(40) 

These expansions are the starting point for justification of the MF theory. 



+ 
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4.2 Restriction of observables 

Let 0( fc ) be an operator acting on the Hilbert space ®j=o^.r Because the creation and 
annihilation operators for bosons are unbounded operators, it is hard to prove that expec- 
tation values of an arbitrary operator calculated by coincides with that by Qmfj- 
Therefore, we now restrict the physical quantities. 

In this chapter, we focus only on the observables 0^ k \a, a*, {/3 r , j3*.}) <G B. The set of 
observables B is defined as follows. All the observables 0^(a, a*, {f3 r , /3*}) € B satisfy the 
following conditions; there exist constants co,c a ,{c r } € C, with liniA_ i , 00 ^ r |c r | 2 < +00, 
{a q > 0}, s > 0, and k > such that 



d m 0( fe )(a,a*,{/3 r ,/3 r *}) 



da ai da a2 ...a 



q=0 



c + c a a + 



(41) 



for all m = 0,1,... and for all {oj}'-L 1 . Here Oj = ±1, a\ := a, and a_i := a*. For 



instance, all the operators expressed by a polynomial of j a, a 1 , {<%} . \ belong to B after 



the correspondence {a/y^N — > a, a) / y/~N — > a*} is made. 

The derivatives d/da and d/da* appear because of the quantum fluctuation, i.e. the 
commutation relation [a, a'] = 1. Roughly speaking, the condition 41 means that we exclude 
observables which are too sensitive to quantum fluctuations. 



4.3 Free time evolution of the 0-subsystem 

Let us define 



—i 



J dt'C p (t') 



U$ : Tex 
For an arbitrary function f(a, a*, {/3 r , /?*}), 

f(a,a*,{pr,fi})U$ = U^f( a (t),a*(t),{(3 r (t),(3* r (t)}), 
where a(t) is the solution of the following equation, 



(42) 



(43) 



d_ 
~dt 

d_ 
V ~dt 



a(t) 



(44) 



Pr(t) = -i (u r /3 r (t) + X r a(t)) , 



with the initial condition a(0) = a and /3 r (0) = j3 r . Because these equations are linear, the 
solutions are written in the form 



a(t) = h a (t) + g aa (t)a + ^2g ar (t)/3 r , 

r 

f3 r (t) = h r (t) +g ra (t)a + ^ g„>(t)(i r >. 



(45) 
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It is noted that 



' \gaa(t)\ 2 +^2\gar(t)\ 2 = l, 

r 

9ra{t)\ 2 + J2\9rr'(t)\ 2 = 1. 



The matrix g(t) is a unitary and symmetric matrix. 

The following property derived from Eq. (43) is important: 



(46) 



u t,s u s,0 — u t,0 u t,s > 



where 



cxp 



-ij*dt>Ci k Xa(t>)) 



(47) 
(48) 



5 Exactness of the MF dynamics in the optical cavity system 
5.1 Statement 

From now on, we prove that the MF theory is exact in the limit of N — > oo in the sense 
that if initially 



JimJ d 2 a (j[ J d 2 ^ Tr li2 ,..., fc O( fc )(a,a*,{/3 r ,/3;}) (q% - Q^ FQ ) = (49) 
for V0( fe ) G 5 and any fixed fe € N, it implies that 

lim^l d 2 a m J d 2 (3 r ) Tr li2 ,..., fc O( fe ) (a, a*, {ft, ^*}) (<* " Qmf,*) = ( 50 ) 



for any fixed time i > 



5.2 Some comments 

We make some remarks here on our result: 

• Our result is quite general: Environments may be attached to the system of interest, 
and the time-dependent driving field may be applied. When an environment is large 
enough, it acts as a thermal reservoir on the system of interest, hence the result of 
this work is applicable to open quantum systems. However, we assumed the special 
Hamiltonian for the O-environment, H-q +H\ in Eq. (21). It is preferable to generalize 
our result to more wide class of O-environmental Hamiltonians in the future. 

• We restricted the class of physical quantities to B in the present work. This restriction, 
however, might be just a technical assumption. To extend to more general class of 
operators is an important issue. Another important issue is to treat the fluctuations. 
For instance, a)a/N = (aa) - 1)/N = (a/ v / iV)(at/ v / ]V) - 1/N belongs to B, but 
(a 1 a — (aJ)(a))/y/~N, whose expectation value converges to a finite value as iV — >■ oo, 
does not. In the seminal work by Hepp and Lieb [16], the authors showed that, 
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under the singular reservoir limit and the approximation to replace the bosons by 
the ensemble of fermions, the equations of motion for the fluctuations are equal to 
the linearized equation of the intensive observables, which is known as the Onsager's 
regression hypothesis [29-31]. We have not been able to confirm the validity of this 
hypothesis without those approximations yet. 

• Without the driving force, and if the thermal reservoirs have an identical temperature, 
the MF dynamics usually predicts thermalization of the system of interest. It means 
that the limit of N — > oo and the limit of t — > oo are interchangeable under such 
an equilibrium situation. In other words, there is no QSS in open systems under 
the equilibrium situation, which is in contrast to closed systems, see [32] for this 
aspect. However, it is less obvious whether or not the thermodynamic limit and the 
long-time limit are interchangeable under the nonequilibrium situation studied in this 
paper. In the present work, we always take the limit of N — > oo first and then the 
limit of t — > oo. In real experimental setup, there can be a situation in which the 
thermodynamic limit ought to be taken after the long-time limit since the number of 
atoms in the cavity is not large enough. 

• We can justify the MF theory, but it is a separate issue whether the MF dynamical 
equations, Eqs. (24) and (25), can be solved. Since the degrees of freedom of environ- 
mental systems remain, it is difficult to solve exactly the MF dynamical equations. It 
is very interesting if we can exactly solve the quantum dynamics of the whole system 
including thermal reservoirs with the help of the MF theory. If it is done, we will 
be able to obtain some insight into effects of the non-Markovian dynamics [25] and 
the accuracy of the perturbative quantum master equation in long times [33, 34] in 
an interacting many body system. When the system of interest is small, which is the 
standard setting of open quantum systems, it is recognized that the non-Markovian 
effect is negligible in the van Hove limit [35]. When the system of interest is large, 
however, it is not obvious whether the use of the Born-Markov quantum master equa- 
tion is justified even for the small coupling between the system of interest and the 
thermal bath because two limits, the thermodynamic limit and the van Hove limit, 
are involved. Therefore it is important to obtain the exact time evolution in a macro- 
scopic open system. It is noted that Merkli and Berman [23] recently obtained the 
rigorous result in this direction for a simple model which is purely dephasing and 
where the relaxation of energy does not take place. 

We hope that the result of this work becomes a good starting point to consider the 
above remaining problems and analyze some interesting nonequilibrium phenomena. 

6 Proof of the statement 

6.1 Finiteness of the average value of observables 

First of all, we prove that all the average values of observables £ B at time t are 

finite as long as they are also finite at initial time t = 0. Unboundedness of operators is 
obstructive for our proof, hence this property is desirable. Once we can show this property, 
the justification of the MF approximation is almost straightforward. 
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From Eq. (39) and Eq. (40), we obtain 
/ d2a (I! / d2 A V.^ 1 " (a,a*,{/3r,/3 r *})Q 



(k) 

N,t 



<B 1 + B 2 + Bz + B A . (51) 



{Bi} are denned by 



x 



B 



Bi = e / d2 « (n / ^] (<*,<**,{&,#» 

JO JO Jo 

/ d ' a (ll J S ^ Tr i,2,-,kO W (a,a*,{/3 r ,/3 r *}) 

x / dtl r^...f^dtjj^^ 



JO 



and 



B * = ^E f d2a (jlf ^ Tr i.2,..,^ (fc) («,«*, {/3 r ,/3 r *}) 

JO JO JO 



1) 



b, = j: 



k + l-l 



i=i 
i-t rti 



N 



J d " a (il / d2 H Tr i,2,..,fcO (fc) (a,a*,{Pr,fi}) 



x [ dt 1 [ 1 dt 2 ...[ 11 ddDgwWC^ 1 ' • • • U^W^Q^ 

v JO JO JO 

Remaining part of this section is dedicated to the evaluation of the upper bounds of Bi 
(* = 1,..., 4). 

Evaluation of Bi 

First we consider the term B\. From Eqs. (34), (35), and (47), we obtain the following 
expression: 



L-l 



1=0 



dti 



E ^ 



(k) — T^i rt k t 1] °,TtMX™ 



{<Tj=±l} 



^—^Tr k+l X k+l U ti:0 Q Nfi 



(52) 
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where xj 1 ^ := Xj and xj 1 - > := xj. By performing integration by part, we obtain 



L-l 



=o 



We put 



i 



X 2^ A fe+l A fc+2 • • • A k+l 

{*j=±l} 



d 



X 







da ci {ti) \da m _ l (ti- 1 ) \ da ai (ti) 







0( k \a(t),a*(t),{f3 r (t),P* r (t)}) 



v r?( fc )lr?( fe + 1 )\ f- r (fe+i-i)\fv(fe+0 / n (fe+0 n (fc+0 
x u t,ti r u ti,t 2 r • • • U t l _ 1 ,t l ( u ti,0 \^N,0 ^MF,0 



(53) 



B 



(k+l) 



t rti 

o Jo 



tl-1 



dtl dt 2 ... Cft, £ < 1 l ) <2-"-</ 

7o K=±i} 



d 



x 



x (a(t), a'(f), {&(*),#(*)}) ^5 K£ j ■ ■ ■ 



fV(fc) 1 ,7(fc+l) \ rV(fc+J-l) \ fr(fe+0 



Then B\ is expressed as 

L-l 



Z=0 



1^ n/^ Tr i. 2 ,.,^ +1) c 



In Eq. (54), the quantity like 



.rarrCO 



appears. We find that 

da 



(k) 

t.s 



da ai (ti) 



da ai {ti)da C2 {t 2 ) ■ ..da an (t n ) 



• , i.,fAk) da ai {t') A (- ai ) x ~( fc) 



0a<n(*i) ^ 



= -ig f dt'U^g aa {t' -1^X^11^ 

Js 3=1 



(As) 
s' 



where X- <T1 ^ X (-) := [xj (•)] is the commutator. Remember that |<7oa(*)| ^ 1 because 



(54) 



(55) 



(56) 



of Eq. (46). Therefore, 



dU, 



(k) 

t.s 



da ai (ti) 



<2kg\\X\\ OQ (t-s). 



(57) 
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Similarly, we can obtain 



(k) 

Us 



da ai (ti)da a2 (t 2 ) ■ ■ ■ da an (t r , 



KpkgWXlUt-sW. 



(58) 



Because of Eq. (58), we can replace dUt,s/da a (t') simply by 2kg\\X\\ 00 (t — s)Ut,s as far as 
it is concerned with the upper bound of the operator norm. This fact makes much easier 
to evaluate the upper bound. 
We define the set 

M,n ■= {{mi,m 2 , ■ ■ ■ ,m n } : 1 < m 1 < m 2 < ■ ■ ■ < m n < I, rrij G N} . 

If one of the elements of _A/i, n is denoted by N^ n = {mi,m 2 , . . . ,m n }, we define Nf n := 
{l,2,...,l}\N hn . 

By using Eq. (58), we then obtain 



B 



(k+l) 



< w ll-X'lloo*)' 



/! 



n 



d 



da„(t) 



0^(a(t),a*(t),{f3 r (t),(3* r (t)}) 



x(k + l- l)(k + l — 2)...(k + n) x (2g\\X\\ 00 t) l - n . 
Because G B, we have 



(59) 



B 



{k+l) 



<^3^ EE E X> 



{aj=±l} n=0 Ni :Tl eMi :n 9=0 



c + c a a(t) + ^2c r p r (t) 



x (k + l — l)(k + 1 — 2) ...(k + n) x (2s||X|| 00 t) , - n . 
Since ^{^=±1} 1 = 2' 5 Z)jv f n eMi „ ^ = ^/( n K^ ~~ n )0; we obtain the upper bound 



(60) 



B 



(k+l) 



< 



E 



1 (fc + Z-1)! 



. — n! (Z - n)!(fc + n - 1)! 



(asllAHoot) 



2l-n 



X E fl ' 
9=0 



c + c a a(t) + ^2c r (3 r (t) 



(61) 



Here we used (k + I — l)!/[(/ — n)\(k + n — 1)1] < 2 k+l ~ 1 and restrict the time interval 
t G [0,*o] where to is determined by 2g>||X||oo*o = 1/2- This restriction will be removed 
later. Then we obtain 



B 



(k+l) 



< 2 k - l - 1 e 2 



K Es 

9=0 



c + c a a(t) + ^2crf3 r (t) 



(62) 
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From Eq. (55) and Eq. (62), we obtain 
B 1 < 2 k e 2K J d 2 a (jj J d 2 ^ 



c + c a a(t) + ^c r /3 r (i) 



Q ( n]o- ( 63 ) 



Here \cq + c a a(t) + Y^ r c r f3 r (t)\ g G B (which will be explicitly shown later). Because at the 
initial time t = 0, for any £ B, the average value is finite by assumption. Therefore, 
Bi is also finite. 

Evaluation of B 2 

B2 is given by 



Bo 



J d2a [PL J d2 ^J V/ 1 {&>#}) J o dh 1^ dt 2 ... 1^ ^ dt L 

u^ k M k J> {k+l) ■ ■ ■ ut^> ik+L - 1] Q { Nt ] 



(64) 



This is written in the form 



where B^ +L ^ is given by 



rt(k+L) l ST^ ( CT1 ) ( CT2 ) v(°l) 

^2 -~9 2^ A fc+l A fc+2 • • • A fc+L 

fo=±l} 

o r a 



dt 2 ... / ^LTr l!2 ,...,fe+L 
Jo 



B 



(k+L) 



n (k+L) 



(65) 







da (TL (t L )\da (TL _ 1 (t L ^i)\ da ai (ti) 
Similarly to the evaluation of 2?i, we obtain the upper bound 



(66) 



B 2 < / d 2 a 



t r t L -i 
dh... cit L L!(2 5 ||X|| 0O ) J 



n/^)i 

x E ^T2' +L - 1 (29iixiut) l - m ^E<' 9 

L 1 

<(2 5 ||X|| 0O t) L £ —2 k+L -\2g\\X\\ O0 t) 



+ c a a(t) + J^c r /3 r (t) 



g=0 
L—m 



(k+L) 
L 



m=0 



X At 



> a max / <i 2 a! 1 I 



d 2 j3 r 



c + c a a(t) + ^2crf3 r (t) 
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Here we again restrict the time interval t G [0, to] (*o was determined by 2g||X|| 00 to 
in the previous subsection). Then we obtain 



B. 



(k+L) 



m=0 



g=0 



co + c a a(t) + ^c r /3 r (t) 



<2- L+k - 1 e 2K ^2< 



q=0 



c + c a a(t) + ^2c r (3 r (t) 



(6 



Substituting it to Eq. (65), we have 



B <2 - L+ k-i y J_ {2K)m y max [ £ a (\\ 

m=0 g=0 \ r 



d 2 /3 r 



Co ~\~ c a a(t) +Y,CrPr{t) 



^^—L+k—l^lK 



> a a max / ci 2 a TT / d 2 /3 r 



o {0) 



c + c a a(t) + ^c r /3 r (t) 



Evaluation of B 3 

Next we evaluate B% that is given by 



5 3 =^E / ^ a (I! / ^) (a,a*,{/3 r ,/3 r *}) 

JO JO JO 



(7' 



If we put 



B 



(k+l) 



{-ig) l t dh dt 2 ... dU 
Jo Jo Jo 

v y(°"i) y(°"2) r («Ti) 

X 2^ A fc+l A fc+2 • • • A k+l 

K=±i} 



x 0<*) (a(t),a*(t), iMt),m) U^}^} • • • U^}*™, (7 



where 



/- k 



for a = l, 



j'=i 



(7 



for<x = -l, 
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we obtain the upper bound of B3 : 
B, 

Similarly to the evaluation of B2, we obtain for < t < to 



h < ^ max I d 2 a (\\ 



d 2 f3 r 



B. 



(k+l) 



W Nt >- 



(73) 



B 



(k+l) 



q=0 



(74) 



Therefore, 



B: 



< ^y {k + l ) 2 k - l ^e 2K ya q max I d 2 a{\\ I d< 



Pr 



c + c a a(t) + y^c r /3 r (t) 



O (0) 

N,ti • 



(75) 



By elementary calculation, we find Ya=\ < 2, and thus 



B 3 < 



c + c a a(t) + ^c r /3 r (t) 



O (0) 
(76) 



Evaluation of S 4 

Finally we evaluate the contribution of H4. Its upper bound is given by 

L 



fi 4 <4y max / d 2 a I TT / d 1 

4 - N^t'e[o,t]J \\1J 



Pr 



B 



(k+l) 



n (k+l) 
Vat// i 



(77) 



where B[ k+l ^ is defined by 



'0 Jo 

X 2^ A fc+l A fc+2 • • • A k+l 

{*j=±l} 



B f +l) : = (-ig)\k + l-l) I dt 1 f tl dt 2 ... 







dti 







X 



d 



da ai {t{) \da (Tl _ l {ti- 1 ) \ da ai {ti) 
x OW (a(t),a*(t), {&(*),#(«)) U&l]^) • • • CS^}" (78) 

We can evaluate the upper bound of B[ k+l ^ by some calculations similar to those in 
previous subsections, and the result is 

q 



B 



(k+l) 



< (fe + /-l)2 fc - i " 1 e 2K ^ 1 



q=0 



co + c a a(t) + ^c r /3 r (i) 



(79) 
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for < t < to- Therefore, we conclude that 
B 4 < ^-2 k ~ 1 (k + l)e 2K Va„ max f d 2 a I TT 



in the time interval < t < to. 



d 2 p r 



c + c a a(t) + ^c r /3 r (t) 



O (0) 
(80) 



Proof of the finiteness of average values of observables 

From the derived upper bounds of Bi (i = 1,2,3,4) in this section, we show that the 
average value of any observable in B at any fixed time t is finite. Collecting the derived 
upper bounds, Eqs. (62), (69), (76), and (80), we obtain 



(0«) t <2 k e 2K ^2a q 

q=0 



Co + c a 



+ 2 k ~ l e 2 * ( 2~ L + 



) > a a max 

) fe[o,t] 



=2 fc e 2K ^, 



q=0 



9=0 



max 

t'e[o,t] 



q=0 

c + c a a(t) + ^2c r /3 r (t) 

r 

c + c a a(t) + ^c r /3 r (t) 



c + c a a(t) + ^2c r (3 r (t) 



(81) 



where b q = 2 k - 1 e 2K (2~ L + (2k + 3)/N)a q . We can see that b q -5- in the limit of N -)• oo 
and L —> oo. Therefore, in the limit of iV — > oo, we can choose b q as an arbitrarily small 

value. Here (-)t ■= f d 2 a(Y[ f d 2 ! /3 r )(-)Qyvt is shorthand notation of the averaged value. In 
order to show the finiteness, we must show that 



max 

t'e[o,t] 



co 



+ c a a{t) + ^c r /3 r (i) 



is finite. 

In order to do so, we put 



and define 



O n {t) :-- 



d + d a a + drfir 



--■■O n , 



d + d a a(t) + d r Pr{t) 



(82) 



(83) 



Because the constants do, d a , and J2 r \dr\ 2 are finite, the generality is not lost if we restrict 
these constants as 

\d\ 2 := \d \ 2 + \d a \ 2 + Y,\dr\ 2 < I- (84) 
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Obviously, O n G B with s = n, cq = do, c a = d a , c r = d r , k = nd a , and a q = 1. 
Now we evaluate the quantity 



X n := max max (O n )f. 
d:\4<it'e[o,t] 



From Eq. (81), we have 



X n < 2 k e 2K (°S)) + 6^ max max (O q (t)) t , , 



(85) 



where b = 2 k - 1 e 2K [2- L + (2k + 3)/iV]. 

Now we use the property of a(t) and /3 r (t). By substituting Eq. (45) into the definition 
of O q (t), we have 



OS) = 



d' + d' a a + d' r /3 r 



(86) 



where 



d'o = do + d a h a (t) + ^2 d r h r (t) =: d- h(t), 

r 

d' a = d a 9aa{t) + ^ d r g ra (t), 
r 

d' r = d a g ar (t) + ^ d r >g r > r (t). 

r> 

The norm of this new vector dl is given by 



\df = K\ 2 + K\ 2 + Y\d' r \ 2 <\h(t)\ 2 + E 



y]d x ig x >S) 



where x, x' G {a, {r}}. Now we use the fact that the matrix g(t) is unitary (absolute value 
of any eigenvalue of g{t) is unity). It yields 

\df<\h(t)\ 2 + l, (87) 

and this is a finite constant for any fixed time t. This also indicates that O q (t) £ £>, there- 
fore, the first term of RHS in Eq. (85) is finite. We write it as f n (t) = 2 k e 2K £™ =0 (0, (t))„. 

We define the vector d" := d'/(\h(t)\ 2 + 1), whose norm is less than or equal to unity, 
|<i 7/ | < 1 regardless of the choice of d. Obviously jd 7/ : \d/'\ < 1 j D \d" : \d\ < 1 j. There- 
fore, from Eqs. (85) and (86), we obtain 

i\ 



max 

e[o,t] 



X n < fS) + b (\h(t)\ 2 + l) V max r 

n 

< fn(t) + b(\h(t)\ 2 + l) V max max 
v y ^d":|d"|<i*'e[o,t] 

n 

= /n(i) + ^(l^)| 2 + l)E^ 



dg + d" a a + £ <A 



(88) 



g=0 
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Since we can choose b as an arbitrarily small number, by starting from Xq = 1, we can 
inductively show that all the X n are finite by using the above inequality. It completes the 
proof of the finiteness of the expectation values, because RHS of Eq. (81) is then finite. 
In particular, we find 



lim lim B2 = lim lim B3 = lim lim B4 = 0. 

L— >-oo iV-^oo L^ooiV->oo L->oo iV->oo 



6.2 Justification of the mean-field approach 



(89) 



The proof proceeds very similarly as the previous section. From Eqs. (39) and (40), we 
obtain the upper bound 



(fc) 

MF,t 



/ d2a (J! / ^,2,.,kO^ (a,a*,{/3 r ,/3 r *}) {q%-Q 

<B[ + B' 2 + B 3 + B 4 



(90) 



Here B3 and B4 are the same ones in the previous subsections. B[ and B' 2 are defined by 



d 2 a 



x 



B':= 



fj y <i 2 /3 r j Th,2,... )fe (fe) (a, a*, {&,#}) 

jT' ^ jT ^ . . . jT -1 ^ W (*)<;> (fe+1) • • • y^+^T ( Q ™ - qSS 

/ Sa (il / ^ Tr i,2,..,fc0 (fc) (a, a*, {&,#}) 
Jo Jo Jo 



r7 .(fc+L-l) w ( fe+ L_i) ( n (k+L) n (k+L)\ 
■■■ U t L -i,t L VV \^N,t L ^MF,t L ) 



(91) 



In the previous section, we showed that limjv->.oo -B3 = lhrijv->oo -B4 = in the time inter- 
val i € [0,to]- Therefore, the proof is completed by showing limjv-^oo B[ = liniAr^oo B' 2 = 
and confirming that we can extend t to all positive values, t £ [0, 00). 
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Evaluation of B[ 

First we consider the term B[. Similarly to the case of B\, we obtain the following expres- 



sion: 



L-l 



B i=E 



1=0 



f^dh dt 2 ... 1 dti Jd 2 a\T[J d 2 p r \ TW,2 ) ...,fc+i 



I 



E v (<Tl) Y (a 2 ) Y (a t ) 
• • • ^k+l 



{<Tj=±l} 

d 



d 



(92) 



It is expressed as 



L-l 



*i = E 



Z=0 



n/^)^A---.^f +,) (o^ ,) -oss5) 



(93) 



where i?| fc+ ^ was defined by Eq. (54). From the assumption of the initial condition (49), 
we can conclude that lirri7v->oo B[ = if we can show B[ k+l) € B. 
In order to show B[ k+l ^ € £>, we must consider the quantity 



d m B\ 



(fe+0 



da T1 da T2 . ..da Tm 



(94) 



where Tj = ±1. 

Calculation analogous to that for B\ yields 



d m B\ 



(k+l) 



da Tl da T2 . ..da Tm 



< 



(sll*lloo*)' 



I 



n 



d 



n 



E E E E E 

{o-j=±l} n=0 N lt „eNi, n P=0 N m ,p£Af m ,p 

d 



da a (t) II J-- 1 - dot T .(t) 



C (fc) (a(t),a*(t)) 



x (fc + l) p (k + l-l)(k + l-2)...(k + n)x {2g\\X\\ 00 t) l ~ n+v . (95) 
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Because 0< fc ) € B, we have 



da T1 da T2 . . . da Tm 



< (g\\X\\ooty 



l\ 



{<7j=±l} n=0 N l} „£Af l}n p=0 N m , P €N m ,p 



x(k + l)P{k + I — l)(k + l — 2)...(k + n)(2g\\X\\ 00 t) l - n+p 



x K n ~ m+P Y,a q 

q=0 



C + C a a(t) + ^ CrPr 



(96) 



Since ^{^=±1} 1 = ^N ln eMi n 1 = ^-/{ n K^ ~ n )D> we obtain the upper bound 



d m B\ 



(k+l) 



da Tl da T2 . . . da Tm 



1 1 



2n-l 



n=0 



[K + 2(k + l)g\\X\\ 00 t] m ^ 



q=0 



co + c a a(t) + ^2c r p r (t) 



(97) 



Jk+l) 



From this expression, we find that B\ £ B for any fixed t. Thus we can conclude that 

(98) 



lim B\ = 0. 



Evaluation of B' 2 

B' 2 is given by 



B' = 



(a,a*,{/3 r ,/? r *}) x J dt^ J dt 2 . . . dt L 



(99) 



This is written in the form 

B> < J Sa (ft J d^)j Tr 1>2 ,..., fc fi^L (q%? + Q<&$) ■ (100) 

where i s given by Eq. (66). From the upper bound oiB { 2 k+L) (see Eq. (68)), it yields 



B> 2 < 2 - L + k -^J d*a (J} / Tr,.., fc £, 



c + c a a(t) + ^2c r /3 r (t) 

r 
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which is almost the same as Eq. (69). Because 



Co + C a a{t) + ^2crP r (t) 



q=0 

is finite and not diverging as L — > oo, it is concluded that 



n (k+L) n (k+L)\ 



lim lim B' 2 = 0, (102) 

L-^oo AT— j-oo 



for < t < to. 



Completion of justification of the mean-field theory 

In the previous subsections, we showed that B\ + B 2 + B% + B4 — >■ as N — > 00 for 
< t < to- It means that Eq. (50) has been derived for < t < to- If we regard t = to as 
a new initial time and repeat the argument, we can show that Eq. (50) is correct for any 
fixed time t > 0. Thus we completes the proof of justification of the MF dynamics. 
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